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Abstract 

This paper presents a general formulation to construct high order numerical schemes by using multi-moment con- 
straint conditions on the flux function reconstruction. The new formulation, so called multi-moment constrained flux 
reconstruction (MMC-FR), distinguishes itself essentially from the flux reconstruction formulation (FR) of Huynh 
(2007) by imposing not only the continuity constraint conditions on the flux function at the cell boundary, but also 
other types constraints which may include those on the spatial derivatives or the point values. This formulation can be 
also interpreted as a blend of Lagrange interpolation Hermite interpolation, which provides a numerical framework to 
accommodate a wider spectrum of high order schemes. Some representative schemes will be presented and evaluated 
through Fourier analysis and numerical tests. 

Keywords: High order scheme, flux reconstruction, multi-moment constraint, nodal formulation, derivative Riemann 
problem, conservation 



1. Introduction 

Solving governing equations point-wisely at the specified nodes (solution points) results in a class of very efficient 
schemes, which are usually known as nodal form, differential form or collocation form, in which no explicit numerical 
quadrature is involved. Examples of this kind are the nodal discontinuous Galerkin (DG) method[3] and the spectral 
collocation methodl^]. In fl, Huynh suggested a general formulation, so called flux reconstruction (FR), from which 
the aforementioned schemes among others can be retrieved from a Lagrange interpolation polynomial with different 
correction functions that ensure the continuity of the cell-wisely constructed flux functions at cell interfaces. The 
DG method can be derived if Radau or Legendre polynomial is chosen for the correction function, while the spectral 
collocation method can be retrieved if the correction function is collocated at the Chebyshev points to a zero function 
inside the cell while matches the modified flux function to the continuous value at the cell boundary. These schemes 
minimize the modification to the primary Lagrange interpolation, and only requires the reconstructed flux function 
itself to be continuous at the cell boundary, which is the necessary condition for the conservation. Wang and Gao 
have implemented the FR in unstructured meshes [ 16], and Vincent et al. have devised some stable subset schemes of 



Huynh's FR formulation with a parameter switching I 15lll4ll . 

If we interpret the flux reconstruction as an interpolation procedure with some constraining conditions given, the 
reconstruction of the modified flux function in Huynh's FR formulation can be essentially viewed as a Lagrange 
interpolation procedure, i.e. all the constraints that lead to the final modified flux function are solely the point values 
(PV) including those at the cell boundaries for conservation and stability and those at the solution points. 

In an alternative direction, we have so far explored the possibility to make use of the Hermite interpolation in 
constructing high order schemes lfl9l \vi\ ITsl I5I l6ll . Compared to the Lagrange interpolation where only PV is used, the 
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Hermite interpolation uses also the spatial derivatives and even integral as the constraints to determine the reconstruc- 
tion interpolation. Hence, we categorize this class as "multi-moment" type schemes against the formulations where 
only one kind "single moment", for example the PV in the Lagrange interpolation, is used. In a recent paper fl, 
we proposed a nodal type formulation, so called multi-moment constrained finite volume method (MCV), where the 
PVs at solution points are updated by evolution equations which are derived via a Hermite interpolation reconstruc- 
tion under the multi-moment constraints including point values, derivatives and integral. The MCV method has been 
implemented for unstructured grids in spherical geometry with competitive performance in numerical accuracy and 
computational efficiency Using multi-moment constraints provides another platform to construct high-order 

schemes. 

We in this paper presents a general formulation that makes use of not only the PV but also other moments like 
derivatives as the constraints to build the modified flux function. The formulation is called multi-moment constrained 
flux reconstruction (MMC-FR), from which new high order schemes can be straightforwardly devised, such as the 
MCV method. 

This paper is organized as follows. The multi-moment constrained flux reconstruction formulation is presented in 
section 2 with comparisons to Huynh's FR formulation. A von Neumann analysis and the convergence rate evaluation 
using grid refinement tests are given in section 3 for some representative schemes of MMC-FR approach, and the 
paper ends with some remarks in section 4. 



2. Multi-moment constrained flux reconstruction formulation 

We consider the following conservation law 

where u is the solution function, and f(u) the flux function. 

The computational domain is divided into / non-overlapping cells or elements Q,- = [jc f _i , jc f+ 1 ], i = 1,2, ■ • • 
and K solution points Xjk,k= 1 , 2, • ■ ■ , K, are set over Q, where the solution k = 1 , 2, • • • , K, is computed. Suppose 
that a proper approximation of the flux function fj(x) is constructed, we can immediately update the solutions within 
Q, by the following point-wise semi-discretized equations at solution points x&c, 



duik 
dt 



dfi(x) 



dx 



1,2, (2) 



The above equation (f2| is of nodal form, also termed as differential form or strong form, which covers a wide range 
of schemes including the collocation methods. The central task left now is how to reconstruct the flux function fi(x). 
In principle, the way to reconstruct fi(x) makes difference among the numerical schemes. 

Given the solution u$ at jCft, k = 1,2, •• • , K, a piece-wise Lagrange interpolation polynomial of degree K—\ for 
Q.i reads 

K 

Ui(x) = 2^ u ik <f> ik (x), (3) 

k=l 

where fakix) is the Lagrange basis function, 



i 1 Tju — Tit 



l=l,M Xik X " 



Given flux f(u) as the function of solution m, we have the flux function at the solution points X&, k — 1 , 2, ■ ■ • , K, 
simply by = f(uik). The piece-wisely reconstructed polynomial for flux function is then obtained as 



fiix) =Y u ffk<Pik{x), (5) 



/>=! 



2 



which is of degree K— 1, same as that for the solution function Uj(x). 

We call (0 the primary reconstruction which is separately constructed over each cell and thus broken from cell 
by cell. The primary reconstruction cannot be directly used to calculate (0, further modification is required to ensure 
at least the C° continuity at the two ends and x i+ i) of cell Q,, which is the necessary condition for numerical 
conservation and computational stability. 

The way to reconstruct the modified flux function fi(x) in an MMC-FR scheme is substantially different from the 
FR of Huynh fl. For comparison and completeness, we describe Huynh's flux reconstruction formulation in 2.1, and 
present the MMC-FR formulation in 2.2 as follows. 

2.1. Huynh's flux reconstruction 

Consider a hyperbolic conservation law, Huynh's flux reconstruction consists of the following steps. 

(i) Given the state variable and the flux function at the solution points, compute the primary reconstruction for the 
flux functions through (0 for all cell Q,, i = 1 , 2, • ■ ■ , /. 

(ii) Compute the flux functions on two sides of the cell boundary x t _i separately by f L , = (;*:,•_ ±) and f R , = 

2 I- J 2 I- J 

(iii) Find the numerical flux function, f® l at the cell boundary by solving 

ff_ k = Riemann^,, (6) 

where "Riemann(- , •)" denotes a solver for the Riemann problem given the values at the left and right sides of 
cell boundary x t _i . 

(iv) The modified flux function over Q, in Huynh's formulation is then computed by 

fax) = /,(*) + [f'f_, - Mx H )]g iL (x) + [f'f +i _ - /i(jt i+ i )]*«(*). (7) 

where gn(x) and giR(x) are the so called correction functions for cell Q, which enforce the continuity at the cell 
boundaries, and satisfies 

gaXx i _ i ) = l; g iL (x iH ) = (8) 

and 

*«(*,_$) = 0; 8iR(x i+ i)=l. (9) 
Both gnix) and giR(x) are polynomials of degree K approximating the zero function, so the modified flux 
function fj(x) is of K degree. Shown in |4], the nodal discontinuous Galerkin (DG) methods [3] and spectral 
collocation method H can be respectively recovered by choosing the correction function to be the Radau 
polynomial and the Lagrange polynomial for the Chebyshev collocation points. See U4J for an comprehensive 
and detailed discussions. 

(v) As long as the modified flux function fj(x) is found, the numerical solutions are updated by the semi-discretized 
equations (fJJ using algorithms for ordinary differential equations, such as the Runge-Kutta scheme. 

2.2. Multi-moment constrained flux reconstruction 

The MMC-FR approach constructs the modified flux function fj(x) in a different manner as follows. 

(i) Compute the primary reconstruction for the flux functions through <j5j f° r a ll ce ll Oj, i = 1, 2, ■ ■ ■ , /. 

(ii) Compute the flux functions on two sides of the cell boundary x t _\ by _f L , = fi-iix^i) and f R , = fi(x t _i), and 

their derivatives by = f™(x,-0 = £ (/-i(^-p) ^d f™? = /^Wp = £ (fii^J)\ 

(iii) Find the flux function, f B , , and its derivatives, / f , at the cell boundary by solving 

I- j XI- 2 

ff_ h = Riemann^,,. ff_^, 
^.DRtemann^,^), 
where "DRiemann(- , ■)" denotes a solver for the derivative Riemann problem (DRP). 
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(iv) The modified flux function f(x) of K degree is then constructed by properly choosing K + 1 constraints of two 
kinds, i.e. a) the continuity conditions of flux function as well as its derivatives at the cell boundaries, referred 
to as boundary constraints, and b) the constraint conditions at some points inside the mesh cell which can be 
directly computed from the primary flux function /;(*), referred to as interior constraints. 

(v) The modified flux function f(x) is then obtained from the K + 1 constraint conditions. Thus, the numerical 
solutions are updated by solving the semi-discretized equations (0 through time integration. 

It is obvious that the major difference between the FR of Huynh and the MMC-FR lies in steps from (ii) to (iv). 
The FR of Huynh only requires the continuity of the flux function at the cell boundaries and all the constraining 
conditions used in determining the modified flux reconstruction are given by the point values, and the interpolation 
is essentially of Lagrange type, whereas the MMC-FR requires the continuities at the cell boundaries of not only the 
flux function itself but also its derivatives, which leads to a Hermite type reconstruction of the modified flux function 
where both the point values and the derivatives are used as the constraints. 

In the MMC-FR formulation, we need compute the derivative Riemann problems at each cell boundary. Prac- 
tically, the approximate Riemann solvers, such as the local Lax-Friedrich (LLF) flux [11] and the Roe's flux ificl 



can be used. The high-order derivative Riemann problems by linear and homogeneous derivative Riemann problems 



are detailed in I113L I12H for the hyperbolic systems. As addressed in 111311 . since the first-instant plays a leading role 
in the interaction of the two states, the derivative Riemann problems with these simplifications provide a reasonable 
accuracy. For the Euler equations in fact, our numerical experiments show that the following linearization to the flux 
functions gives adequate accuracy in terms of both numerical error and convergence rate, 

f = Au, j\ - An..-- - . fl m] =Au™, (11) 

where A is the Jacobian matrix obtained by A = df/du. Provided the derivatives of the state variable u from the 
cell-wise reconstructions for both sides of a discontinuity, one can find the derivative flux of any order at the expense 
of the conventional Riemann problem. 

In the present study, the reconstruction is carried out in terms of flux function / itself, instead of the conservative 
variable u. For example, the numerical approximations for flux function / and its derivatives at cell boundary x lp can 
be calculated by 



J xip <-) W xip J xip l P x f 1 ip x xip xip J I 

f (12) 



where we make use of the relation: u = A~ l f = (RAR l y l f = RA l R l f with A being the diagonalized matrix of 
the eigen values. The eigen matrices R and R l , as well as the eigen values in A are directly evaluated by the point 
values at Xi P . 

2.3. Implementation of the MMC-FR formulation 

As discussed above, the MMC-FR uses not only the point value of flux function but also its spatial derivatives 
as the constraints, and thus is more flexible with greater freedom to experiment with. We in this section presents the 
general procedure to construct the modified flux with some concrete examples. 

Given K solution points, the primary reconstruction fix) is of K - 1 degrees. The modified flux reconstruction 
f(x) should be at least of K degree to get the solutions of (f2]i u± of degree K — 1. So, we need K + 1 or more 
constraints for the reconstruction. Suppose we use the continuity constraint conditions up to the Ar^th derivative of 
the flux function at the left boundary and the k^th derivative of the flux function at the right boundary, we have the 
boundary constraints as 

4r 1 (* i -p=/I ml f. m = 0,l,--,k L ; 

1 (13) 

Am\, x Am]B n 1 ; 

f xi • m = 0,l,--- ,k R ; 

where the flux function is included as the case of Oth order derivative, i.e. f(x) = 
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We leave the rest constraints, kj — K - ki - + 1 in number, to be the interior constraints and determined by 
coinciding the modified flux function with the primary flux function /;(jc) at some points (referred to as constraint 
points) in terms of point value or spatial derivatives, i.e. 

/!r'W) = W). fc' = l, ■••,&/; m' = 0,l,--. (14) 

Note that the constraint points x^ are not necessarily the same as the solution points X&. 

From the constraint conditions given in (TPTb and (fT4l . the K degree modified flux function fi(x), is uniquely 
determined. 

Next, we give some examples of deriving the flux functions by the MMC-FR approach. For brevity, we use a local 
coordinate system ^ e [-1, 1] that transforms the real mesh cell x e [Xi-1/2, ^1+1/2] by 

X — Xj_l 

f = 2— r — -1. (15) 

Ax; 

where Ax ; = x Mjl - Xi-1/2. 

The constraints cTT~3T > and (fT4l are recast correspondingly as 



(16) 



and 

Am 



/f ] (D = ff S (D = (|) ; J™, m = 0,l,-,k R ; 

f^'Xfa) = f\f\fa), k' = \,---,ku m'=0, (17) 



where fa = 

The time evolution equations (f2]l for updating the solutions are correspondingly 

dj^ = _idj\idjm\ , = 1)2> ...^. (18) 

• Three-point scheme 

We consider a scheme having three solution points, fi, £2 and £3, where the solution u&, k = 1,2,3, are 
updated every time step. We assume that the flux f(u) is a function of solution u. The values of the flux at the 
corresponding points, fy, k - 1, 2, 3, are computed directly. The primary flux function is then built by 

3 

/«(£) = £ /***(£)• (19) 
k=\ 

where 

3 

n t~t (2o) 

is the basis function of the Lagrange interpolation. 

The continuous flux function and its first order derivative at the cell boundaries in (TToT l are computed from steps 
(i)-(iii)of2.2. 

The modified flux function is constructed by using boundary constraints up to ki — k^ — 1, while no interior 
constraint is used, kj = 0. That is, the modified flux function of degree 3 is determined by the following 
constraint conditions, 
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ik-i) = /f(-D; 
/;(!) = /f(i); 
4, 1] (-D = 4, 1]S (-D; 

L/1 ] (D=4 1]S ( 1 )- 



(21) 



(22) 



(|2TT > is a Hermite interpolation to determine the modified flux function which is written in a polynomial form 

as, 

+^ (3/f (i) - 3/f (-i) - /|(-i) - 

+~ (2/f (1) + 2f'f(-\) + /|(-1) - /|(1)) . 

The first order derivative (gradient) of (OTl reads then, 

3 



m) =4 (/?(-d - /f tt) + /|(-d + /|d))^ 2 

+2 (3/; s (i) - 3/f(-D - /|(-d - /|(d) . 



(23) 



The derivatives of the modified flux function at the solution points are obtained as 



( djm \ 

I # la 



/f£2 - /#(&); 



/f/3 - /fi(ft)- 



(24) 



The solutions are then immediately computed by ( fTST l with a proper time integration algorithm. 
Choosing different solution points results in different schemes. We give two examples next. 

- Equidistant points 

If the collocation points are equally spaced and including the cell ends, i.e. £i = -1, £2=0 and £3 = 1, we 
retrieve the third-order MCV scheme^, where the left/right-most solution points coincide with the cell 
boundaries. In this case, the continuity conditions of flux function at the cell boundaries are automatically 
satisfied, and only the derivatives of the flux function need to be computed from the DRP. 
The derivatives of the modified flux function at the solution points are obtained as 



fm =/^i) = /|(-i); 

fa = ~ (3/f (1) - 3/f (-1) - /|(-1) - /|(1)) ; 

JfB= /fife) = /|(D. 



(25) 
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It is noted that in the three-point case equidistant points are identical to the Chebyshev-Gauss-Lobatto 
points. It is straightforward to show the following conservation property, 



i: 



f (g)d€ = J] (wkffft) = \f m + \ka + \f f o = jf (D ~ jf (-1), 



(26) 



where = J_j <pk(^)d^ are the weights of numerical quadrature. 
- Chebyshev-Gauss points 

We use the Chebyshev-Gauss points, i.e. ft = - V3/2, ft = and ft = y3/2, as the solution points. In 
this case, the continuity conditions of flux function and its derivatives are not necessarily satisfied. So, all 
of them have to be computed from the DRP. 

The derivatives of the modified flux function at the solution points are obtained as 



fyi = ^ (3/f(D - 3/f (-1) + (5 + 4 V3)/|(-l) + (5 - 4 V3)/|(l)) ; 
fea = \ (3/f(D - 3/f(-l) - /|(-1) - /|(1)); 
/f/3 = ^ (3/f (1) - 3/f (-1) + (5 - 4 V3)/|(-l) + (5 + 4 V3)/|(l)) . 



(27) 
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Again the numerical conservation yields from 
>i 3 



(28) 



k=\ 



Four-point scheme 

When we use four solution points, ft, ft, ft and ft, the primary flux function is constructed by the Lagrange 
interpolation, 



fffi = 2^(ft- 



(29) 



fc=i 



After continuous flux function and its first order derivative at the cell boundaries in ( fT~6b are computed, the 
modified flux function is constructed by using boundary constraints up to = = 1, while one interior 
constraint is used, kj = 1 . The modified flux function of degree 4 is determined by the followings, 

{ ik-D = /f(-i); 



^•(i)=y; s (i); 
4 1] (- 1 ) = /| 1]S (- 1 ); 
4 1] (!) = 4 11B(1); 



(30) 



U(0)=/i(0). 

(f30b is a hybrid expression including both Hermite (the first four) and Lagrange (the last one) interpolations, 



and can be written in a polynomial form as, 

7;(ft =2 (-2/f (-D - 2/f(D + 4/(0) - /|(-1) + f%Q))? 

e3 



4 

+ \ (4/?(-l) + 4/f (1) - 8/,(0) + /|(-1) - /|(1))^ 2 

+ 7 (3/f (i) - 3/f (-i) - /|(-1) - + ./X0). 



(31) 



The first order derivative (gradient) of (|3~T1 reads then, 

'hm = (-2/; s (-i) - 2/; s (d + 4/ko) - /|(-D + f'§{i))f 
+1 (/f (-d - ffd) + /f (-i) + /|(i))^ 2 

+ 5 (4/f (-1) + 4/f (1) - 8/KO) + /|(-1) - /|(l))^ 
+ 7 (3jf (i) - 3/fC-l) - /|(-D - /|(1)) . 



(32) 



The derivatives of the modified flux function at the solution points are obtained as 

\ la 
l dfi(g) \ 
\ ) B 

U k 

The solutions are then immediately computed by (TT~ST > with a proper time integration algorithm. 
An alternative to the constraint conditions (l30b is 

M-i) = /f(-D; 
/*(!)= if (i); 





= &(fi); 


/fi2 


= /«»(&); 


/f«3 


= .4(6); 


/fi4 


= 



(33) 



/J ] d) = 4/ ls (i); 
cfiMO) _ d 2 fm _ /[2](0) 



(34) 



Here, we retain the curvature of the primary interpolation rather than the value at the cell center. The modified 
reconstruction for the flux function is then, 

+\ (/fa) - /f (-D + 4, 1]S (-i) + % v <x))t + ~/JW 
+ \ (4/f (D + 4/f (-d + /|(-i) - /|(-d - 4 21 (0)) . 



(35) 



The first order derivative (gradient) of (1351 at the solution points can be directly computed in the same way. We 
call this scheme the MCV4 with central 2nd derivative constraint, MCV4_C2D in short. 

The formulations discussed above work for different collocation points. We give some examples as follows. 

■1, £2 = -1/3, £3 = 1/3 and £4 = 1, which results in the 



- Equidistant point 

We use 4 point equally spaced points, £1 



4-order MCV scheme[6], where the left/right-most solution points coincide with the cell boundaries. In 

8 



this case, the continuity conditions of flux function at the cell boundaries are automatically satisfied, and 
only the derivatives of the flux function need to be computed from the DRP. 

The derivatives of the modified flux function of the MCV4 scheme ( 1321 at the solution points are obtained 

as 



[fax =#(-D; 

~ka = ^ (2/f(D - 34/f (-1) - 8/|(-l) - /|(1) + 32/f (0)) ; 
/fo = ^ (34/f (1) - 2./f(-l) - 8/|(l) - /|(-1) - 32/f (0)) ; 



(36) 



27 

l/f«4=/|(l). 

For the MCV4_C2D scheme, the modified flux function ( f35l > results in different formula, 



7«n 


= /|(-i); 








= ^(i8/f(i)~ 


i8/f(-i) 


- 4/|(-l) - 5/|(l) - 8/| ]S (0)) ; 




= l(i8./f(i)- 


i8/f(-i) 


-4/|(l)-5/|(-l) + 8/| ]S (0)); 




= /|(D- 







(37) 



We can directly prove the following numerical conservation for both (l36"T l and (l37l >. 



fc=i ^ 



1 « 



3 ~ 



3 ~ 



1 s 



4 -/ f ii + j/ja + + 4/^4 = /f(l) - ^(-1). 



(38) 



- Chebyshev-Gauss-Lobatto points 

For Chebyshev-Gauss-Lobatto points, %\ = -1, £2 = -1/2, £3 = 1/2 and £4 = 1, the derivatives of the 
modified flux function ( 1321 ) at the solution points are obtained as 

fpi=f§(-l)l 

ka = (~3/f (1) " 21/f (-1) - 3/|(-l) + /|(1) + 18/f (0)) ; 
^ (21/f (1) + 3/f (-1) - 3/|(l) + /|(-1) - 18./f (0)); 



(39) 



/f<3 
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l/f/4 = /|(l). 

Those for the MCV4_C2D scheme can be obtained from d35l > directly. Again, the following relation of 
conservation holds for MCV4 and MCV4_C2D, 



k=\ v 



1 * 



1 s 



(40) 



- Chebyshev-Gauss points 

If one chooses Chebyshev-Gauss points, the four solution points locate inside the cell element £1 = 
cos(77r/8), £2 = cos(57r/8), £3 = cos(37r/8) and £4 = cos(7r/8), the derivatives of the modified flux function 
at the solution points can be obtained immediately by (13T1 >. For reader's convenience, we write them in a 
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decimal form as 

' f m = - 0.160763093/f(l) - 0.380433007/f(-l) + 0.0635243020/|(l) + 0.7168057838/|(-1) 
+ 0.541 196100/f(0); 

fga = ~ 0.0131 164397/f(l) - 1.293446526/f (-1) - 0.0048660179/|(1) - 0.2754640679/|(-l) 
+ 1. 306562965/^(0); 

fga =1.293446526/f(l) + 0.0131 164397/f(-l) - 0.2754640679/|(l) - 0.0048660179/|(-1) 
- 1.306562965/f(0); 

f m =0.380433007/f(l) + 0.160763093/f(-l) + 0.7168057838/|(1) + 0.0635243020/|(-l) 



• 0.541 196100/f(0). 



(41) 



The numerical conservation is verified by the following equality with round off error, 

4 ' \ 

fai&dt;) = 0.2642977396/ f! -i + 0.7357022609/ f/2 + 0.7357022607/ fi3 + 0.2642977393/ f/4 



k=l x 



= /?(!)" /?(-!)• 



Semidiscrete Spectra 
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(42) 



Figure I: The spectrum of the semi-discrete schemes. 



Five-point scheme 

When we use five solution points, ft, ft, ft, ft and ft, the primary flux function is constructed by the Lagrange 
interpolation, 



(43) 



k=\ 
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In this case, we are given more freedom in choosing the constraints to construct the modified flux function. We 
show two examples as follows. 



,5. All constraints are imposed 



- MCV5 scheme 

Equidistant solution points are located at & = - 1 + 2(k — 1)/4, k = 1,2, 
on the cell boundary in terms of the flux derivatives up to ki-k^- 2, 

f/X-i) = ./f(-i); 

Ml) = f(l); 

4- ] (-i) = 4- ]S (-D; 

4/ ] (D=4; ]S (i); 

/? ] (-D = 4 2] Vi); 

l4 2] (D=4? ]S (i). 

(l44l i is a Hermite interpolation which leads to the next polynomial, 

' fco = ^ (-3/f(-D + 3/f (i) - 34 11s (-d - 34, 11B (i) - 4, 2]S (-d + 4 2is (i))^ 5 

+^(4/ ls (-i)-4 1]S (i)+4f ]S (-i)+4 2]S (D)^ 

+ 1 (5jf (-1) - 5/f (1) + 54^(-l) + 54/^(1) + 4^(-l) - jf\l))f 

+ \ H4/ ] Vd + - 4 2]s (-d - 4 2]B (i))^ 2 

+ ^ (-15/f (-D + 15/f (1) - 74, 1]S (-1) - 74, 1]S (1) - 4f ]S (-l) + 4f ]S (l))£ 

+?g (s/f (-D + s/f (i) + 54, 1]S (-d - 54, 1]S (i) + 4 2]S (-d + 4 2]S (d) . 

The derivatives of the modified flux function at the solution points then read 

fei = 4/ ls (-i); 

^ (135/f (1) - 135/f (-1) + 8l4 Ils (-l) - 95.4 1]B (1) + 274, 21S (-1) + 2l4, 2is (l)) ; 
^ (15.^(1) - 15./f(-D - 74, 1]S (-1) - l.§ YB (l) - 4f ]S (-l) + 4 2]S (D) ; 
-i- (l35/f(l) - 135/f(-l) - 954, 1|S (-1) + 8l4; ]S (l) - 2l4, 21B (-l) - 274, 21S (1)) ; 



fan. = 

/f/3 = 

fm - 



256 

[fsa = 4 11S (D. 



The numerical conservation can be obtained from the following relation 

5 / M 



£ - ^ (7/jfl + 32/^-2 + + 32/ fi4 + 7/^5) = f'f(l) - f'f(-l). 



(44) 



(45) 



(46) 



(47) 



- MCV5_PV24 scheme 

We locate the solution at the Chebyshev-Gauss-Lobatto points, & = - cos((k - 1)tt/4, k = 1, 2, • • • ,5, 
and combine the multi-moment constraints at the cell boundary in terms of the flux derivatives up to 
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Table 1: Numerical errors and convergence rate for three-point schemes. 



Scheme 


(0 = 7r/8 


oj — 7r/16 


order 


MCV3 


-3.25 x 1(T 4 - 3.33 x KT 5 / 


-2.06 x 10- 5 -1.07x lQ- b i 


2.99 



kh — k-R — 1 and two collocation constraints at point k — 2 and k — 4 inside the cell, 

f/k-D = /f(-D; 
/KD = /f(i); 
4- ] (-i) = 4- ]S (-i); 
4 1] (i) = 4/ ]S (D; 
/(&) = /<(&); 
J i (Z4) = f i m. 

(|4B1 is a mixture of Hermite interpolation and point collocation (Lagrange interpolation) from which the 
modified flux function is obtained in a polynomial form as, 



(48) 



m = 



5 (/f c-i) - fm) + \ (4 1]S (-d + 4- ]S (i)) + 2 V2 + /-(f 4 )) 
2 (-/f(-D - /fa)) - \ (-4/ ] Vd + 4, 1]S ( 1 )) + 2 (/»(&) + /«(&)) t 



2 



4 (/f c-d + /fa)) + ^ (4 11S (-!) - 4' 18(1) ) - 4 ^(&) + 
J (/f (-D - /fa)) + \ (4/ ls (- 1 ) + 4 1IS(1) ) + 2 ^ (-/■(&) + /■(&)) 
- § (/f (-d + /fa)) + \ - 4" s (- 1 )) + 2 V5(-/i(6) + / { (f4)) ■ 

The derivatives of the modified flux function at the solution points then read 

'fen =4- ]S (-l); 

/„, = (| - 2 VSjjf (1) - (I + 2 Vlj/f (-1) + I ( V2 - 1)4^(1) -I(V5+ i)4, I]S (-D + 



V2 



(7/K6) + /*(&)); 



/fa =5 U S (-!) " /fa)) + J (4/ ]S (D + 4, 1IS (-D) + 2 V5 (/■(&) - /i(£0) ; 

fm = (1 + 2 V5) /f (i) - g - 2 V5) jf (-1) + 1 ( V2 - i)/^(-i) -i(V2+ i)4/^(i)- 



V2 

— (7/;(&) +/;(&)); 

>=$ ]S a). 

The numerical conservation can be obtained from the following relation 

5 



2 £ = (/f« + 8 -fe + 4 /f« + 8/f« + 7 fe) = /f (1) - /f (-D- 



(49) 



(50) 



(51) 
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- MCV5_2D24 scheme 

Another choice is to use the 2nd-order derivatives of the primary interpolation function as the constraints 
at points k = 2 and k — 4, 

' J(-i) = /f(-i); 
Mi) = ./f(D; 



4 ! 11 (- 1 ) = 4« 1!S (-i); 



4/ ] (i)=^ra); 

df 

d 2 m) d 2 fm 



(52) 



fg ] m 



de de / ™ ( ^ 4) " 

Constraint conditions d52l results in another scheme, so called MCV5_ 2D24. In the same manner, the 
derivatives at the solution points, fy^, k = 1, • • • ,5, can be directly obtained , and the numerical conserva- 
tion can be proven. 

Shown above, the MMC-FR formulation can be implemented in an efficient manner that only involves the nodal 
values. Higher order schemes can be straightforwardly devised by increasing the solution points and proper constraint 
conditions. As we will show later, although convergence rate of at least Kth order can be easily achieved for a 
Appoint scheme, the numerical errors and the stable CFL restriction differ by the constraint conditions. Moreover, 
the numerical conservation can be rigorously guaranteed as long as the modified flux function at cell boundaries are 
continuous. 



Table 2: Same as Tabled but for four-point schemes. Note that the errors are examined at a) = jr/4 and a) = 



Scheme 


CD 


= 7r/4 


CO 


= 7T/8 


order 


MCV4 


-4.91 x 10" 


s + 9.42 x 10- s z' 


-7.88 x 10" 


'+3.17x10-^ 


4.02 


MCV4.C2D 


-1.95x10" 


4 + 5.71 x 10" 4 j 


-3.15 x 10" 


6 + 1.91 x \Q~H 


3.96 



3. Fourier analysis 

In this section, we evaluated the numerical schemes previously discussed by examining the Fourier mode trans- 
ported with the following advection equation, 

du du „ ,^„ N 
¥ + (53) 
The theoretical tools used hereafter are mainly developed in J2,3l. 
We use a wave solution, 

u(x,t) = e IbKx+t) , (54) 

and represent it on a grid whose ith cell is defined by [i— 1/2,/+ 1/2] where the solution points x/j, k — 1 , 2, • • • , K, 
are located. The solutions are then u& = e Ia>( -'^^ 2 \ We consider a scheme constructed over cell / with the information 
from its upwinding cell (/ - 1). Recall that boundary flux and its derivatives, f' B and /j'" 1 " 8 , are completely upwinding 
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and / = u in this particular case, the schemes discussed above are summarized by 



I dun \ 

dt 
dun 

dt 

du iK 
dt 



C(i-l)12 
C(/-l)21 C (; _i)22 

V C(f-i)Kl C(i-\)K2 



CQ-\)\K Cm C,12 
C(i-l)2K C/21 C/22 

C(i-\)KK CiKl CiK2 



Ci\K ^ 
Ci2K 

CiKK ) 



From the wave solution, we have = e /w ('+ft/ 2 ) and Un-m = e Ia> Uik, which reduces 



W(j-1)2 
U(i-\)K 
Uq 

V UiK ) 

into 



(55) 



dug 

dt 
du-a 

dt 

du iK 
v dt 



( c ( /_i) U e lb> +cai c ( ,-i)i 2 e lt0 + c a2 
C(,-i)2ie _/ " + Cq\ c (! -_ I)2 2e _/ " + c i22 



C(i-\)K\e l0J + c iK i c (i -i )K2 e luJ + c iK2 



CQ-\)\Ke Il ° + CaK \( «il ^ 



CQ-i) 2 Ke Iw + CaK 



C(i-\)KKC + C,KK ) 



\ u iK J 



(56) 



or 



dUi 

-di = SU - 



(57) 



As discussed in 0], the constraint conditions make the essential difference in matrix S, while different arrange- 
ments of the solution points result in similar matrices which have the same eigenvalues. In the following discussions, 
equi-distanced points are used for the MMC-FR schemes. 



Numerical dispersion 
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Figure 2: Numerical dispersion relations for the schemes. 



The properties of the numerical schemes can be examined by analyzing the eigenvalues of ( l57b . FigQ]shows the 
spectrum (collection of all eigen values) of S for different schemes. It is observed that all eigenvalues lie on the left 
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half of the real axis, i.e the negative real parts indicate that all the schemes are stable under the CFL conditions. The 
allowable CFL numbers for computational stability can be estimated by the largest eigenvalue, the spectral radius p for 
each scheme, i.e. a scheme has a larger spectral radius has to use a smaller CFL number for computational stability. 
We know from Fourier analysis that Pmcv3 = 6.0, Pmcva - 9.78 and Pmcva.cid - 5.42. The MCV4_C2D scheme has 
a spectral radius even smaller than the three-point MC V3 scheme. Shown later, our numerical tests for pure advection 
equation verify the observations from the spectral radius analysis presented here. 




Figure 3: The numerical dissipation relations for the schemes. 



The numerical errors of different schemes can be examined by comparing the principal eigenvalue of S, with 
the exact solution, -Ico, of the advection equation (l53l for initial condition, 

u(x,0) = e Iwx . (58) 

The error of a given semi-discrete formulation is 

E(co) = + Ico, (59) 

and the convergence rate is evaluated by 

m = \\n\ ,5 ,L I / ln(2)| - 1 ■ (60) 



E(lo/2) 

The numerical errors of the MCV3 scheme are given in Table [1] The MCV3 scheme has a 3rd-order convergence 
rate. The errors of the four-point schemes are given in Table [2] The two MMC-FR schemes, MCV4 and MCV4_C2D, 
are stable and of 4th-order convergence rate. 
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The dispersion and dissipation relations of the spatial discretizations can be evaluated by plotting the real and 
imaginary parts of the principal eigenvalues as a functions of the wave number a>. From Figf2] we find that all schemes 
agree well with the exact solution for a> < 4, which reveals a superior numerical dispersion compared to conventional 
finite difference method or finite volume method. For the four-point schemes, the MCV4 scheme adequately recovers 
the dispersion relation for all waves of u> < 7. The MCV4_C2D scheme, however, is less accurate for high wave 
numbers. Fig|3] shows that both MCV4 and MCV4_C2D are more accurate in numerical dissipation even for high 
wave number. 

We evaluated the convergence rates of the schemes discussed above by solving the linear scalar equation (l53l l on 
gradually refined grids with the smooth initial condition defined by 

u(x, t = 0) = sin(7rx). (61) 

A periodic boundary conditions are specified. The normalized errors ft, (2 and („ at t — 2 are given in Table [3] 

All schemes get the expected convergence rates. The numerical errors, however, varies according to the constraint 
conditions imposed. It is observed that the constraints on the point values at the interior points will reduce the 
numerical errors, whereas the constraints in terms of the second order derivatives result in larger numerical errors. 



Table 3: Normalized errors and convergence rates of the MCV type schemes. 



MCV3 


Mesh 


ft 


1 1 -order 


ft 


^2-order 






„ -order 




10 


1.29 x 10- 2 




1.43 x 10- 2 




2.03 x 10- 2 








20 


1.69 x 10~ 3 


2.93 


1.87 x 10~ 3 


2.93 


2.69 x 10- 3 




2.92 




40 


2.14 x 10~ 4 


2.98 


2.37 x 10- 4 


2.98 


3.37 x 10- 4 




3.00 




80 


2.68 x 10~ 5 


3.00 


2.98 x 10~ 5 


2.99 


4.22 x 10~ 5 




3.00 


MCV4 


Mesh 


ft 


1 1 -order 


ft 


^2-order 


too 


(\ 


„ -order 




10 


2.06 x 10~ 4 




2.26 x 10~ 4 




3.24 x 10~ 4 








20 


1.31 x 10~ 5 


3.98 


1.46 x 10" 5 


3.95 


2.09 x 10~ 5 




3.95 




40 


8.32 x 10- 7 


3.98 


9.25 x 10- 7 


3.98 


1.31 x 10- 6 




4.00 




80 


5.24 x 10~ 8 


3.99 


5.82 x 10~ 8 


3.99 


8.25 x 10~ 8 




3.99 


MCV4.C2D 


Mesh 


ft 


1 1 -order 


ft 


£2 -order 




L 


„ -order 




10 


1.22 x 10~ 3 




1.33 x 10~ 3 




1.91 x 10~ 3 








20 


7.84 x 10~ 5 


3.96 


8.75 x lfT 5 


3.93 


1.26 x 10" 4 




3.92 




40 


5.00 x 10- 6 


3.97 


5.56 x 10- 6 


3.98 


7.90 x 10~ 6 




4.00 




80 


3.15 x 10~ 7 


3.99 


3.50 x 10~ 7 


3.99 


4.96 x 10~ 7 




3.99 


MCV5 


Mesh 


ft 


1 1 -order 


ft 


f^-order 


too 


(\ 


„ -order 




10 


5.21 x 10~ 5 




5.72 x 10- 4 




8.18 x 10 -5 








20 


1.67 x 10- 6 


4.96 


1.85 x 10~ 6 


4.95 


2.65 x 10~ 6 




4.95 




40 


5.28 x 10~ 8 


4.98 


5.86 x 10 s 


4.98 


8.31 x 10~ 8 




5.00 




80 


1.65 x 10-" 


5.00 


1.84 x 10-" 


4.99 


2.61 x 10-' 




4.99 


MCV5.2D24 


Mesh 


£1 


t\ -order 


ft 


£ 2 -order 


t m 




„ -order 




10 


4.58 x 10~ 5 




5.02 x 10- 4 




7.20 x 10~ 5 








20 


1.46 x 10- 6 


4.97 


1.62 x 10- 6 


4.95 


2.32 x 10- 6 




4.96 




40 


4.63 x 10~ 8 


4.99 


5.13 x 10- 8 


4.98 


7.29 x 10- 8 




4.99 




80 


1.46 x 10-" 


4.99 


1.62 x 10-' 


4.98 


2.29 x 10-' 




4.99 


MCV5.PV24 


Mesh 


£1 


t\ -order 


ft 


{ 2 -order 


t m 




„ -order 




10 


3.48 x 10- 6 




3.85 x 10~ 6 




5.48 x 10~ 6 








20 


1.07 x 10- 7 


5.02 


1.18 x 10- 7 


5.03 


1.70 x 10~ 7 




5.01 




40 


3.33 x 10-" 


5.01 


3.70 x 10-' 


5.00 


5.25 x 10-' 




5.02 




80 


1.06 x 10-'° 


4.97 


1.17 x 10-'° 


4.98 


1.66 x lO^ 10 




4.98 



The largest allowable CFL numbers for different schemes are also evaluated through the numerical tests for ad- 
vection equation with a third order Runge-Kutta time integration scheme. The results are shown in Table|4] It is found 
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that using the point values of the primary interpolation at the interior points tends to reduce the maximum stable 
CFL number, while using the 2nd-order derivatives can increase the stable CFL number. Consistent with the obser- 
vation from FigJT] which show a smaller spectral radius of MCV4_C2D compared to MCV3, the 4th-order scheme 
MCV4_C2D can use larger CFL number than the 3rd-order scheme MCV3. 



Table 4: Largest allowable CFL numbers. 



Schemes 


MCV3 


MCV4 


MCV4X2D 


MCV5 


MCV5.2D24 


MCV5.PV24 


CFL(max) 


0.425 


0.275 


0.485 


0.21 


0.235 


0.165 



4. Concluding remarks 

We have presented a general formulation using the multi-moment constraints to construct high-order schemes for 
hyperbolic conservation laws. The formulation accommodates a wide range of schemes under the unified framework 
of flux reconstruction. 

Leaving the practical schemes to be further explored for specified applications, we give some analysis in this paper 
to show the basic feature of the multi-moment constrained schemes, and obtain the following observations: 

• The multi-moment constrained finite volume schemes presented in |0] are stable and possess numerical accuracy 
superior to the conventional finite volume method. 

• The numerical conservation can be rigorously guaranteed as long as the modified flux function at a cell boundary 
is shared by the two neighboring cells. 

• The location of the solution points is not sensitive to the numerical results. 

• Constraints of the point values at interior points improve numerical accuracy, but tends to suffer a more restric- 
tive CFL condition for computational stability. 

• Constraints in terms of the 2nd-order derivatives, i.e the curvature of the primary reconstruction, greatly relieve 
CFL restriction for computational stability. It is possible to construct schemes of higher order convergence rate 
that allow larger stable CFL number. 
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